Under consideration for publication in J. Plasma Phys. 



1 



Computationally efficient methods for 
modelling laser wakefield acceleration in the 

blowout regime 

B. M. COWAN 1 !, S. Y. KALMYKOV 2 }, A. BECK 3 
X. DAVOINE 4 , K. BUNKERS 2 , A. F. LIFSCHITZ 5 , 
E. LEFEBVRE 4 , D. L. BRUHWILER 1 , 
B. A. SHADWICK 2 , and D. P. UMSTADTER 2 

^ech-X Corporation, 5621 Arapahoe Ave. Ste. A, Boulder, CO 80303, USA 

2 Department of Physics and Astronomy, University of Nebraska - Lincoln, NE 68588-0299, 

USA 

3 Centrum voor Plasma- Astrofysica, Departement Wiskunde, Katholieke Universiteit Leuven, 
Celestijnenlaan 200B, B-3001 Leuven, Belgium 

4 CEA, DAM, DIF, Arpajon F-91297, France 

5 Laboratoire d'Optique Appliquee, ENSTA, Ecole Polytechnique, CNRS, 91761 Palaiseau, 

France 

(Received ?; revised ?; accepted ?. - To be entered by editorial office) 

Electron self-injection and acceleration until dephasing in the blowout regime is studied 
for a set of initial conditions typical of recent experiments with 100 terawatt-class lasers. 
Two different approaches to computationally efficient, fully explicit, three-dimensional 



particle-in-cell modelling are examined. First, the Cartesian code VORPAL (Nieter & Cary 



2004 ) using a perfect-dispersion electromagnetic solver precisely describes the laser pulse 
and bubble dynamics, taking advantage of coarser resolution in the propagation direc- 
tion, with a proportionally larger time step. Using third-order splines for macroparticles 
helps suppress the sampling noise while keeping the usage of computational resources 
modest. The second way to reduce the simulation load is using reduced-geometry codes. 



In our case, the quasi-cylindrical code c alder- CIRC ( Lifschitz et a/.||2009| ) uses decom 



position of fields and currents into a set of poloidal modes, while the macroparticles move 
in the Cartesian 3D space. Cylindrically symmetry of the interaction allow using just two 
modes, reducing the computational load to roughly that of a planar Cartesian simulation 
while preserving the 3D nature of the interaction. This significant economy of resources 
allows using fine resolution in the direction of propagation and a small time step, mak- 
ing numerical dispersion vanishingly small, together with a large number of particles per 
cell, enabling good particle statistics. Quantitative agreement of the two simulations indi- 
cates that they are free of numerical artefacts. Both approaches thus retrieve physically 
correct evolution of the plasma bubble, recovering the intrinsic connection of electron 
self-injection to the nonlinear optical evolution of the driver. 



1. Introduction 

Relativistic Langmuir waves driven by short, intense laser pulses in rarefied plasmas 
maintain accelerating gradients several orders of magnitude higher than those accessible 
in conventional metallic structures ( Tajima fc Dawson|1979 Gorbunov fc Kirsanov|1987 
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Esarey et a/.|2009 ). The technical simplicity and compactness of these laser-plasma accel- 
erators (LPAs) is attractive for a broad range of applications, such as nuclear activation 
and on-site isotope production (Leemans et al. 2001 Reed et al. 2007), long-distance 
probing of defects in shielded structures (|Ramanathan et a/.||2010[ ), and testing radia- 
tion resistivity of electronic components (Hid ding et a/.||201l| ). Realisation of compact, 
inexpensive, bright x- and gamma-ray sources using electron beams from LPAs (Rousse 



et a/.1|2004[ |2007| |Ta Phuoc eroT|2005| |Kneip et a/.|2010| |Cipiccia et aJ]2Qll\ holds the 
promise to enable a much wider user community than can be served by existing large- 
scale facilities. These applications are not especially demanding as regards electron beam 
quality, and in fact sometimes draw benefits from poor beam collimation and a broad en- 
ergy spectrum (Hiddin g et a/.|2011 ). However, there are also important applications with 
much tighter beam requirements. Such applications include generating coherent x-rays 
using an external magnetic undulator (iGruner et a/.||2007| |Schlenvoigt et a/.||2008fr|a 



Fuchs et a/.||2009|) , producing x-rays for the phase contrast imaging (Fourmaux et al. 
2011 Kneip et a/.||2011 ), and high-brightness, quasi-monochromatic gamma-ray Comp- 



ton sources ( Leemans et a/.|[2005 Hartemann et a/.||2007 ); these require electron beams 



with a mufti- kA current, low phase space volume, and energy in the few-gigaelectronvolt 
(GeV) range. 

Achieving this high level of accelerator performance is a major near-term goal of the 
LPA community. Modern laser systems capable of concentrating up to 10 Joules of energy 
in a sub-50 femtosecond pulse ( Yanovsky et a/.|2008 Froula et a/.|2009 Kneip et a/.|2009| 



Fourmaux et al. 2011[) make it possible to achieve the so-called blowout (or "bubble") 



regime, which is desirable due to its technical simplicity and scalability (Gordienko 
Pukhov|[2005l |Lu et aZ.||2QQ7| ). In this regime, motion of the electrons in the focus of the 



laser pulse is highly relativistic. The laser ponderomotive force expels plasma electrons 
from the region of the pulse, while the fully stripped ions remain essentially immobile, 
creating a column of positive charge in the laser wake. The charge separation force 
attracts bulk plasma electrons to the axis, creating a closed bubble devoid of electrons. 



This co-propagating electron density bubble (|Rosenzweig et a/.||1991 Mora & Antonsen 


1996 Pukhov & Meyer-ter-Vehn|2002 


Gordienko & Pukhov|2005 


Lu et al. 2006 


guides 


the laser pulse over many Rayleigh lengths ( 


Mora & Antonsen 


1996 Lu et al. 


2007 


)• 



The bubble readily traps initially quiescent background electrons, accelerating them to 
hundreds of megaelectronvolts (MeV) over a few mm, creating a collimated electron 
bunch (Pukhov & Meyer-ter-Vehn 2002 Kalmykov et al. 2011a). It is in this regime 



that the first quasi-monoenergetic electrons were produced from laser plasmas in the 



laboratory (Geddes et al. 2004 Mangles et al. 2004 


Faure et al. 


2004), and the GeV 


energy range was approached (Leemans et a/.||2006| Karsch et a/.||2007 Hafz et al. 


2008 


Froula et al. [2009; 


Kneip et a/.||2009| |Clayton et a/.||2010| Lu et a/.||2011t Liu et al. 


2011 


Pollock et a/.||2011 





Multi-dimensional particle-in-cell (PIC) simulations have played a key role in under- 
standing the physics of the fully kinetic, strongly relativistic blowout regime. The PIC 
method ( Hockney fc Eastwood||198Tj Birdsall fc Langdon||1985 ) self-consistently models 



both electromagnetic fields and charged particles, representing field quantities on a grid 
and particles in a continuous phase space. Given sufficient computing power, electromag- 
netic PIC codes can simulate the plasma electrons (and ions, if necessary), the laser pulse 
driving the plasma wake, and the dynamics of electrons injected into the accelerating po- 
tential. In particular, two- and three-dimensional PIC simulations have been essential 



in understanding the dynamical nature of the electron self-injection process (Xu et al. 
20051 IQguchi et a/.|2008irWu et a/.|20091|Zhidkov et al. |20 101 [Kalmykov et a/.|2009H20101 
2011a|6|c ). However, to capture precisely the correlation between driver dynamics, elec- 
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tron self-injection, and GeV-scale acceleration in the bubble regime, a simulation must 
meet a number of challenging requirements. 

Optimisation of a GeV-scale LPA performance, even with the use of massively par- 
allel computation, is a challenging task especially because of the necessary cm-scale 
laser-plasma interaction length. The laser energy is used most effectively if electrons are 
accelerated until they outrun the bubble and exit the accelerating phase, at which point 
they will have gained the maximum possible energy in an LPA stage, 

^«2.7 7 ^ 3 P^MeV. (1.1) 
Acceleration to this dephasing limit occurs over the distance ( Lu et a/T]|20Q7 ) 



L d « 0.6A 7 s 8/3 j Ptw- (1-2) 

Here, Ptw is the laser power in terawatts (1TW = 10 12 W), j g = uJo/uj pe ^> 1 is the 
Lorentz factor associated with the linear group velocity of the pulse in plasma, c^o is 
the laser frequency, Ao = 2ttc/ujo is the laser wavelength, uj pe = (noe 2 /m e eo) 1 / 2 is the 
electron Langmuir frequency, m e is the electron rest mass, no is the background electron 
density, e is the electron charge, and eo is the permittivity of free space. The scalings 
and (1.2) imply that the pulse remains self-guided , namely, that it remains longer 



than cu;^ 1 ( Sprangle et a/.|l990 Gorbunov et a/.|[2QQ5 ), and that its power exceeds the 
critical power lor relativistic self- focusing, F CY = I6.27 2 GW (|Sun et a/.|l987). Increasing 

oe 



the electron energy therefore requires reduction of the electron plasma density, increasing 
both the bubble velocity and size, 

£ acc « 0.9A o7 ' /3 Ptw. (1-3) 
where L acc is the length of the accelerating phase of the wakefield (roughly equal to the 

bubble radius). Electron dephasing scales as ~ 4 ^ 3 and thus the final energy gain 

2/3 

scales as Ed ~ n . For instance, reaching 1 GeV energy with a 200 TW pulse and 
a wavelength of Ao = 0.8 um may be achieved in a 0.47 cm length plasma with density 
no = 3.5 x 10 18 cm -3 , and doubling that energy would require nearly four times the 
plasma length and three times lower density, also increasing the bubble size by ^40%. 
Simulations of LPA commonly use a moving- widow, where the simulation box propagates 
with the speed of light colinearly with the laser pulse. This optimisation notwithstanding, 
even the experiments with currently operating 100 TW systems bring forth the task of 
modelling the pulse propagation in cm-length plasmas, with the size of 3D simulation 
box on the order of hundred (s) of microns longitudinally and transversely. 

The greatest challenge arises from the great disparity of physical scales between the 
laser wavelength and plasma length, which is the hallmark of high-energy laser-plasma 
acceleration. The need to resolve the laser wavelength, Ao ~ 1 pm, fixes the grid resolu- 
tion, and, due to stability conditions ( Courant et a/.|[l967 ), also limits the time step to 



a small fraction of uj^ 1 . Furthermore, the strong localisation of the injection process im- 
poses even stricter limit on grid resolution; the vast majority of injection candidates are 
concentrated in the inner lining of the bubble (the sheath), and penetrate into the bubble 
near its rear, where the sheath is longitudinally compressed to a few tens of nanometres 
( Wu et a/.|2009 Kalmykov et a/.|2011a ). Resolving this structure, together with ensuring 



sufficient particle statistics in the sheath, is necessary to avoid excessive sampling noise 
and eliminate unphysical effects. In this situation, extending the plasma length to cen- 
timeters and increasing the size of the simulation window to hundreds of microns, while 
at the same time maintaining sufficient macroparticle statistics, would require solving 
Maxwell's equations on meshes amounting to billions of grid points, and advancing 1- 
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10 billion macroparticles over millions of time steps. Performing such simulations with 
standard electromagnetic solvers and particle movers requires a national-scale supercom- 
puting facility. As a result, an attempt to reproduce the long time scale evolution of the 
laser and the bubble together with fine details of the electron self-injection dynamics 
is usually a compromise between affordable simulation load and unavoidable coarseness 
of the results. However, the high precision of modern LPA experiments and high beam 
quality requirements of the applications are rather unforgiving to these compromises and 
do not tolerate numerical artefacts ( |Cormier-Michel et a/.| [2QQ8). 

These considerations make it clear that PIC algorithms must be modified in order 
to reduce the required computational resources without compromising precision. One of 
the main directions is development of electromagnetic solvers that minimize numerical 
error while using the lowest possible grid resolution. One particular limitation of PIC 
that requires high longitudinal resolution is that of numerical dispersion. In PIC, elec- 
tromagnetic fields are typically updated using the finite-difference time-domain (FDTD) 
method on a staggered Yee grid (Yee 1966 Taflove & Hagness 2QQ5[ ). This method is 
second-order accurate, and since it is explicit and local, it parallelizes efficiently enabling 
large-scale simulations. However, it is known that this algorithm experiences numerical 
dispersion error for waves propagating along the axis, which leads to errors in the group 
velocity of the laser pulse. This artificial slowdown of the driver and the bubble leads 
to incorrect dephasing of accelerated electrons and also permits synchronisation of the 
sheath electrons with the bubble, leading to their unphysical injection. Mitigating this 
effect by using higher resolution increases the computation time quadratically. Because 
of the deleterious effects of numerical dispersion in FDTD schemes, efforts have been 
made to develop perfect dispersion algorithms, which exhibit no numerical dispersion for 
waves propagating along a grid axis. For accelerator applications, several modifications to 
FDTD have been described that correct for numerical dispersion using implicit methods 
( Zagorodnov et a/.] 2003; Za gorodnov fc Weiland|200"5| ) . Because LPA simulations tend to 
be quite large-scale (using thousands of processor cores), an explicit algorithm is desir- 
able for reasons of computational efficiency. Such an algorithm has been described in 2D 
( Pukhov[T 999 ) and in 3D for cubic cells (Karkka inen et a/.|2006 ). These algorithms have 
also been explored for LPA as a means of reducing noise in boosted-frame simulations 
( Vay et al. ||20lT| ). 

In this paper, we use two complementary simulation codes (with different numerical 
approaches and physics content) to explore physical phenomena involved in self-injection 
and acceleration of electrons until dephasing under typical conditions of recent experi- 
ments with 100 TW-class lasers. We use a newly-developed perfect-dispersion algorithm 
(Cowan et al. 2012) implemented in the fully explicit 3D Cartesian VORPAL simulation 
framework ( Nieter fc Cary|2004 ) , subsequently referred to as vorpal-pd. The algorithm, 
briefly described in Sec. [2j eliminates numerical dispersion in the direction of pulse prop- 
agation. Thus, even with a relatively large longitudinal grid spacing (~15 grid points per 
Ao), the correct group velocity of a broad-bandwidth laser pulse is obtained. 

The other code used here, calder-CIRC, uses cylindrical geometry. This code uses 
poloidal mode decomposition of fields and currents defined on a radial grid, while macro- 
particles retain their full 3D dynamics in Cartesian coordinates ( |Lifschitz et al. 2009[ ). 
Well-preserved cylindrical symmetry of the laser-plasma interaction enables using just 
a few lower-order modes. Neglecting higher-order, non-axisymmetric contributions to 
the wakefields and currents makes it possible to approach the performance of a 2D code. 
CALDER-CIRC thus allows for fast, extra-high resolution runs with excellent macroparticle 
statistics ( |Kalmykov et al. ||2010[ |2011a[ ) 



The paper is organized as follows. In section [2] we outline the main features of the 
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recently implemented perfect dispersion algorithm in the vorpal-pd code. Section [3] is 
dedicated to the benchmarking of vorpal-pd against calder-CIRC. Sec. [4] summarizes 
the results and indicates the directions of future work. 



2. The perfect dispersion method 

In this section we give a brief overview of the perfect dispersion method we use; a more 



complete description together with detailed benchmarks will be presented in ( Cowan et al. 
2012). Our method is based on that in ( Pukhov|[T9 99 ; Karkkai nen et a/.||2006 ), in which 
the FDTD algorithm is modified by smoothing the fields in the curl operator in one 
of Maxwell's equations. We choose to smooth the electric fields for the magnetic field 
update; our update equations are then 



D t B = -V' X E, D t E = c 2 V XB 



J 



(2.1) 



where J is the electric current deposited from particle motion. Here D t is the finite 
difference time derivative, Vx is the standard finite difference curl operator, and V'x 
is the modified curl operator. Our modification to the curl operator involves applying 
smoothing transverse to the coordinate axis along which the derivative is taken. For 
instance, when computing dE y /dx, E y is smoothed in the y and z directions. This is 
equivalent to applying a smoothing operator before the numerical derivative operator. 
The electric field is smoothed only for the update of the magnetic field; the smoothed 
fields are not stored for the next time step. 

The smoothed curl operator V' X is formed by modifying the finite difference operation. 
If Di is the numerical derivative operator in the i-th direction, then for the modified curl 
we use DiSi in place of Di, where Si is the smoothing operator for the derivative. The 
smoothing operator S x is defined by the stencil in the y and z directions 



lyz Pz 
Py OL x 
lyz Pz 



lyz 
lyz. 



(2.2) 



and similar relations hold for cyclic permutations of the coordinate indices. The coeffi- 
cients c^, fa, and jij are chosen to guarantee that waves propagating along the x axis 
(the laser propagation direction in our simulations) in vacuum experience no numerical 
dispersion, as described in ( Cowan et a/.||2Q12 ). The only constraint is that the longitu- 
dinal grid spacing Ax must satisfy Ax < Ay, Az for the transverse grid spacings Ay and 
Az. 



3. Benchmarking 

While a technological path to high-quality GeV beams exists, experimental progress is 
impeded by an incomplete understanding of the intrinsic relation between electron self- 
injection and nonlinear optical evolution of the driver, and hence by the lack of suitable 
criteria for selection of the optimal regimes that produce beams with the smallest possible 
phase-space volume. Control and optimisation of the fully kinetic, intrinsically 3D process 
of electron self-injection is a daunting task. It involves a systematic study of the links 
among the dynamics of self-injection and the nonlinear optical processes involving the 
laser pulse and the bubble. 

Due to the extended acceleration length, the interaction of the laser pulse with the 
plasma is rich in nonlinear phenomena. Even a Gaussian beam which is perfectly matched 
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to the electron density gradient in which it propagates is not immune to nonlinear op- 



tical processes. Oscillations of the pulse spot-size due to non-linear refraction (Oguchi 



et al. 2008; Zhidkov et al. 2010 Kalmykov et al. 2010), self-phase modulation leading 



to the formation of a relativistically intense optical piston (Tsung et a/.||2002| |Lontano 



fc Murusidze] [20031 |Faure et at] [20051 |Pai et a/.|[20T0l |Vieira et a/.||2010| |Kalmykov 
et a/.||2011a|6| ), and relativistic filamentation ( Andreev et a/.||2d07 Thomas et a/.||2007[ 



2009| ) are processes which result in pulse deformations. Electron self-injection appears to 
be extremely sensitive to such changes in pulse shape, which lead to contamination of 



the electron beam with polychromatic, poorly collimated background (Kalmykov et al. 



20116). Such contamination is readily seen even in simulations with idealised initial con- 



ditions (Kneip et a/.||2009 Froula et a/.||2009| |Martins et a/.||2010 Kalmykov et a/.||2010 



2011a). The complicated modal structure of the incident pulse further aggravates the 
situation, leading to continuous off-axis injection, collective betatron oscillations ( |Glinec| 
et a /.|200 8 ; Mang les" et a/.|2009 Cummings fc Thomas|20lTl ) , and electron beam steering 



(Popp et al. 2010). In practice, these phenomena currently preclude operation reliable 



enough to enable high-precision user experiments; reported islands of stability for self- 
injection in laser and plasma parameter space remain relative ly narrow (|Karsch et al. 



20071 [ Mangles et a/.|2007| |Thomas et al. |2007| |Hafz et a/. (20081 |Maksimchuk et a/.|2008] 



Wiggins a/.|20T0 ). Numerical codes used in predictive modelling of LPAs must be able 
to reproduce these phenomena with high precision in order not to confuse the instability 
of acceleration caused by physical processes with unphysical artefacts caused by intrinsic 
deficiencies of numerical algorithms, such as numerical dispersion, high sampling noise, 
and grid heating. 

3.1. Simulation parameters 
The simulations presented here extend the earlier case study by ( Kal mykov et a/.|2011a 



and use the same set of initial conditions. A transform-limited Gaussian laser pulse with 
full width at half-maximum (FWHM) in intensity tl = 30 fs, wavelength Ao = 0.805 um, 
and 70 TW power is focussed at the plasma border (x = 0) into a spot size ro = 13.6 um, 
and propagates in the positive x direction. The laser pulse is polarised in the y direction. 
The peak intensity at the focus is 2.3 x 10 19 W/cm 2 , giving a normalised vector potential 
of ao = 3.27. The plasma density has a 0.5 mm linear entrance ramp followed by a 
2 mm plateau and a 0.5 mm linear exit ramp. The density in the plateau region, no = 
6.5 x 10 18 cm -3 , corresponds to j g ~ P/P C r ~ 16.3 and dephasing length w 1.7mm. 

The simulations carried out with vorpal-pd use grid spacings of Ax = O.O6A0 = 
48.3 nm longitudinally and Ay = Az = 0.5A = 403 nm transversely, with four macropar- 
ticles per cell. Use of third-order splines for the macroparticle shapes reduces the sampling 
noise, mitigating the adverse effect of the coarse grid. The domain in the vorpal-pd 
simulation is 72 um long and 91 um wide, and is surrounded transversely by a 16-layer 
perfectly-matched layer absorbing boundary. The code is fully parallelised, and was run 
using 6 144 cores on the Hopper supercomputer at the National Energy Research Sci- 
entific Computing Center (NERSC). Completion of a typical run took ~3 x 10 5 CPU 
hours. 

The c alder- CIRC simulation uses 45 macroparticles per cylindrical cell, formed by 
the revolution of the grid cell around the propagation axis. The longitudinal grid spacing 
is Ax = 0.125c/co>o ~ 16 nm. The aspect ratio Ar/Ax = 15.6 (where r = \Jy 2 + z 2 ) : 
and the time step At = 0.1244CJQ -1 . With these grid parameters, numerical dispersion 
is negligible, and sampling noise is significantly reduced. This high resolution simula- 
tion does not indicate any new physical effects compared to the vorpal-pd simulation, 
and does not exhibit significant differences in the quantitative results. Well preserved 
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cylindrical symmetry during the interaction (confirmed in the VORPAL-PD simulation) 
enables us to approximate fields and currents using just the two lowest-order poloidal 
modes, thus reducing the 3D problem to an essentially 2D one. These results confirm 
the earlier established fact (Lifschitz et al. 2009) that, in the case of a linearly polarised 



laser, higher-order modes contribute only weakly to the electric field. Comparison with 
the results of the VORPAL-PD runs shows that our restriction to only two modes is suffi- 
ciently precise to reproduce all relevant physical effects, and to simulate the propagation 
through a 3 mm plasma in 2 625 CPU hours on 250 cores. 

3.2. Formation of quasi-monoenergetic bunches and physical origin of dark current 

Upon entering the plasma, the strongly overcritical pulse rapidly self-focuses, reaching 
its highest intensity at x ~ 0.8 mm, soon after entering the density plateau. Full blowout 
is maintained over the entire propagation distance. In both simulations, electrons are 
accelerated until dephasing in two distinct stages, each characterised by completely dif- 
ferent laser pulse dynamics. Transverse evolution of the laser pulse is the hallmark of 
Stage I. The pulse spot size oscillates, first causing expansion and then contraction of 
the bubble. The bubble expansion produces self-injection of electrons from the sheath; 
stabilisation and contraction of the bubble extinguish injection, limiting the beam charge 
to a fraction of a nC. Phase space rotation creates a well-collimated quasi-monoenergetic 
bunch long before dephasing. Further acceleration (Stage II) is dominated by longitudinal 
(temporal) self-compression of the pulse, leading to gradual elongation of the bubble and 
continuous injection, producing a polychromatic, poorly collimated energy tail with a 



few nC charge. This two-stage evolution has been noticed in earlier simul ations (| Froula 
et a/.|[2009 Kneip et al. 2009), and explained in detail in ( Kalmykov et a/.||2011 



The correlation between the plasma bubble evolution and the self-injection process 
is quantified in figure [I] Panel (a) shows the length of the accelerating phase on axis, 
viz. the length of the region inside the bubble where the longitudinal electric field is 
negative. Panel (b) shows the longitudinal "collection phase space", viz. momenta of 
macroparticles reaching the dephasing point, p x (x = ^deph), vs. their initial position in 
plasma. Panel (c) shows the collection volume: the initial positions of electrons reaching 
the dephasing point. Comparison of these three panels shows that electrons are injected 
only during the periods of bubble expansion. 

During Stage I, radial oscillation of the laser pulse tail inside the bubble causes al- 
ternating expansion and contraction of the first bucket, clearly seen in the progression 
from x = 0.6 to 1.24 mm in figur e p^ a). The bubble size oscillates around the average 
value predicted by the estimate ( |1.3| ), L acc & 9.5 um. Electron self-injection into the 
oscillating bubble leads to the formation of a quasi-monoenergetic component in the en- 
ergy spectrum. At the end of Stage I, at z « 1.24 mm, the bubble contracts to the same 
size in both runs, truncating the tail of injected bunch and expelling electrons injected 
between x = 0.825 and x ~ 0.95 mm. These electrons do not reach dephasing and thus 
are missing in figures [TJb) andjljc). Electrons injected between x = 0.65 and 0.825 mm, 
remain in the bubble and are further accelerated. This well-separated group of particles 
is clearly seen in figure [ljb) . In both the VORPAL-PD and calder-CIRC simulations, 
these electrons reach dephasing first, preserving low energy spread, and are accelerated 
to the highest energy, E w 500 MeV. The bubble expands more rapidly and stabilises 
sooner in the VORPAL-PD simulation, causing stronger reduction of the phase velocity 
in the subsequent buckets (second and third). Hence, in contrast to the CALDER-CIRC 
run, VORPAL-PD gives a noticeable amount of charge trapped and preaccelerated in these 
buckets. These electrons, indicated by the red ellipse in figure [TJb), are swallowed by the 
expanding first bucket during Stage II and are further accelerated, contributing to the 
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Figure 1. (a) Length of the accelerating phase vs. propagation distance in CALDER-CIRC (black) 
and vorpal-pd simulations (red/grey). Expansion and contraction of the bubble due to non- 
linear focusing of the driver (Stage I) is followed by continuous expansion caused by pulse 
self- compression (Stage II). (b) Longitudinal momentum of electrons reaching the dephasing 
point, Xdeph ~ 2.4 mm, vs. their initial longitudinal positions. Black dots are the CALDER- 
CIRC macroparticles; the colourmap represents the normalised number density of vorpal-pd 
macroparticles. Electrons are injected only during periods of bubble expansion. A quasi-mo- 
noenergetic bunch forms during Stage I and maintains its low energy spread until dephasing, 
indicated by the group of early-injected particles with E « 500 MeV. Groups of electrons en- 
compassed by the ellipses were injected into the second and third buckets, to be further captured 
and accelerated by the expanding first bucket. Continuous injection during Stage II creates a 
polychromatic energy tail, (c) Collection volume: initial radial offsets of electrons reaching de- 
phasing limit Rm = v yf n + zf n , vs. their initial longitudinal positions x- lYl . Black (red/grey) dots 
are CALDER-CIRC (vorpal-pd) macroparticles. This collection volume indicates that the vast 
majority of electrons are collected from a hollow conical cylinder with a radius slightly smaller 
than the local bubble size. 



dark current. This contribution, however, appears to be fairly minimal in comparison to 
the amount of continuously injected charge during Stage II. 

The leading edge of the laser pulse constantly experiences a negative gradient of the 
nonlinear index of refraction. As a result, by the end of Stage I, it accumulates con- 
siderable redshift. During Stage II, plasma-induced group velocity dispersion slows the 
red-shifted spectral components relative to the unshifted components, leading to the front 
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Figure 2. Electron density (in cm -3 ) in the plane of laser polarisation in CALDER-CIRC (top 
row) and vorpal-pd simulations (bottom). Panels (a) and (d) show the fully expanded bubble 
in the middle of Stage I, (b) and (e) the fully contracted bubble at the end of Stage I, and (c) 
and (f) the bubble in the vicinity of electron dephasing point at the end of Stage II. x = ct is 
the trajectory of the laser pulse maximum in vacuum; (a) and (d) correspond to the distance 
x = ct ~ 930 um from the plasma edge, (b) and (e) to x = ct ~ 1210 um, and (c) and (f) to 
x = ct « 2364 um. Before the dephasing point, the bubble, elongated and deformed due to the 
laser pulse self-compression, traps considerable charge. Beam loading, however, is yet unable to 
terminate self-injection (cf. panels (c) and (f)). 



etching and pulse self-compressing into a relativistically intense, few-cycle long optical 



piston ( Tsung et a/.|2002 Lontano fc Murusidze|2003 Faure et a/.|2005 Kalmykov et al 
2011a| ). As the pulse transforms into a piston, the bubble constantly elongates, resulting 



in copious trapping and creating a poorly collimated, polychromatic tail, clearly seen 
in figure []Jb). At the dephasing point, ^deph ~ 2.4 mm, the bubble size becomes nearly 
twice the estimate L acc w 9.5pm based on the scaling law (1.3). Even though figure[IJa) 
shows a larger bubble expansion in the vorpal-pd run, the sections of collection phase 
space corresponding to Stage II look nearly identical for both codes in figure [TJb) . 

The collection volume depicted in figure [TJc) indicates that the electrons are collected 
from a conical shell with a radius slightly smaller than the bubble radius. This structure 
of the collection volume indicates that the vast majority of trapped and accelerated 
electrons have impact parameters of sheath electrons ( Tsung et al. [2006 Wu et a/T]|2009 



Pukhov et a/.||2010||Kalmykov et a/.||2010||2011a| ). Collection volumes in the vorpal-pd 



and calder-CIRC runs are almost identical during Stage I, whereas the radius of the 
cone is larger for vorpal-pd during Stage II, on account of the greater expansion due 
to pulse diffraction. 

Snapshots of electron density, longitudinal phase space, and energy spectra at the 
points of maximal expansion and contraction of the bubble are presented in figures |2j 
[3j and [4] Data for panels (a), (b), and (c) in these figures are from the calder-CIRC 
simulation, and for panels (d), (e), and (f) from the vorpal-pd simulation. 

The fully expanded bubble in the middle of Stage I is shown in figures [2^a) and [2^d) . 
As soon as the bubble fully expands, injection terminates. Uninterrupted injection of 
sheath electrons before this point produces a large spread of longitudinal momentum 
and energy, shown in figures [3^a), |3|d), and[4^a). 
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Figure 3. Electron density (arbitrary units) in longitudinal phase space in CALDER-CIRC (top 
row) and vorpal-pd simulations (bottom). Each panel corresponds to the same panel of fig- 
ure [2] Full expansion of the bubble saturates injection and initiates phase space rotation (panels 

(a) and (d)). Contraction of the bubble terminates injection, clipping the rear of injected bunch, 
eliminating low-energy tail. Phase space rotation makes the bunch quasi-monoenergetic (panels 

(b) and (e)). Elongation and deformation of the bubble due to the laser pulse self-compres- 
sion causes continuous injection, producing an electron beam with a continuous spectrum of 
longitudinal momenta (panels (c) and (f)). 
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Figure 4. Electron energy spectra in c alder- circ (black) and vorpal-pd simulations 
(red/grey). Panels (a), (b), and (c) correspond to the phase space snapshots (a) and (d), (b) 
and (e), and (c) and (f) of figure p] respectively, (a) At the point of full expansion, the electron 
energy spectrum is broad, (b) Full contraction of the bubble suppresses the low-energy tail and 
reduces the energy spread. Electrons from the second bucket contribute to the background, seen 
in the diffuse peaks around 150 MeV. (c) Continuous injection caused by the bubble expansion 
and deformation produces a massive polychromatic tail. The leading bunch, at E « 500 MeV, 
reaches dephasing, but is still distinct from the tail. 



The slight contraction of the bubble between x = 0.95 and 1.24 mm truncates the 
bunch. Electrons injected at the very end of the expansion interval are expelled, while 
particles remaining in the bucket are further accelerated. The transverse self-fields of 
the bunch are unable to prevent the bucket contraction. Snapshots of the contracted 
bubble are presented in figures [2^b) and[2|e). During the contraction interval, the tail 
of electron bunch, exposed to the highest accelerating gradient, equalises in energy with 
earlier injected electrons, thus producing a characteristic 'IT shape in the longitudinal 
phase space. This feature (also observed in the similar situation by (Lu et al. 2007)) is 
clearly seen in figures [3jb) and|3|e). As a result of this evolution, quasi-monoenergetic 
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^5mono E mono Ai/ mono ^N,y &N,z 

calder-circ 214 215 20 6.87 7.08 
vorpal-pd 193 245 20 10.7 6.08 

Table 1 . Parameters of the quasi- monoenergetic bunch (E > 200 MeV) at the end of Stage I 
(cf. the spectra in figure H[b)). Qmono is the charge in pC; Emono is the energy corresponding to 
the spectral peak (in Mev); AE mono is the absolute energy spread (FWHM) in MeV; £n, v and 
£n,z are the normalised transverse emittance (in mm mrad) in and out of the laser polarisation 
plane, respectively. 



bunches form in both VORPAL-PD and CALDER-CIRC simulations at the end of Stage I. 
These quasi-monoenergetic spikes with < 10% energy spread can be seen in figure Qb). 
In addition to the quasi-monoenergetic spikes, these energy spectra also reveal diffuse 
features near 150 MeV, corresponding to the electrons trapped in the second bucket; these 
particles can be seen in the snapshots of electron density shown in figures [5|b) and[2|e). 
These electrons, however, never equalise in energy with the leading high-energy bunch. 

Parameters of the bunches, summarised in table [I] appear to be very similar. Nor- 
malised transverse emittances presented in this table are calculated according to the 
usual definition e Nii = (m e c)- 1 [((pf ) - (Pi) 2 )((r 2 ) - (r^) 2 ) - {{pin) - {ri){pi)) 2 ) 1/2 , where 
i = y and z correspond to the emittance in and out of polarisation plane. The beam 
asymmetry is more pronounced in the VORPAL-PD simulation, presumably on account 
of the inclusion of the complete electromagnetic field, in contrast to just two poloidal 
modes in CALDER-CIRC. 

Agreement between the two codes worsens during Stage II. As has already been noted, 
the bubble expansion is larger in the vorpal-pd simulation. As a result, the amount of 
continuously injected charge at the dephasing point (2.5 nC) is about 60% higher and 
the divergence of the continuously injected beam (80 mrad) is about twice that in the 
CALDER-CIRC simulation. The difference in charge can be easily inferred from figure [4jc). 
On the other hand, parameters of the leading bunches are in reasonable agreement, with 
the central energy 485 ± 20 MeV in vorpal-pd against 515 ± 25 MeV in calder-CIRC 
run. In both simulations, the emittance of the quasi- monoenergic component increases by 
~30% over its value at the end of Stage I. The lower energy of the leading bunch in the 
VORPAL-PD run can be explained by its earlier dephasing due to more rapid expansion 
of the bubble. 

Both codes agree that the bubble not only elongates during Stage II, but becomes 
more and more asymmetric in the laser polarisation plane. The "pennant-like" bubble 
shape is responsible for massive off- axis injection, leading to the noticeable beam centroid 
oscillations in the laser polarisation plane seen in figures |2|c) and|2|f). Such phenomenon 
has been observed in similar situations by others ( Glinec et a/.||2008| ). This violation of 
symmetry is a manifestation of carrier-envelope phase effects in the interaction of a 



relativistically intense, linearly polarised, few-cycle piston with the plasma (Nerush & 



Kostyukov 2009). Conversely, in the plane orthogonal to the laser polarisation, both 



the bucket and the beam remain perfectly symmetric (not shown). Surprisingly, the two 
poloidal radiation modes of CALDER-CIRC still capture the field evolution well. Inclusion 
of higher order modes should improve the situation. On the other hand, figures |2jc) and 
[2jf ) indicate that electromagnetic solvers of both codes agree on the group velocity of 
the laser pulse even in the situation where the pulse shrinks down to less than two cycles 
and remains strongly relativistic. This means that (a) poloidal mode decomposition does 
not damage dispersion in the axial direction, and (b) the coarse grid and dispersion 
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properties of vorpal-pd are sufficient to describe well the extreme case of pulse spectral 
broadening to Auj ~ ujq and compression to nearly a single cycle. 

Examination of the bubble evolution and collection volumes (cf. figure [TJ, together with 
individual snapshots of electron density in coordinate and longitudinal phase space, indi- 
cate that, in spite of the great difference in the algorithms, vorpal-pd and calder-circ 
reproduce the same correlation between the evolution of the bubble and the self-injection 
of sheath electrons, and agree quantitatively on the parameters of quasi-monoenergetic 
beams produced by the oscillating bubble. Self-injection begins, terminates, and resumes 
at exactly the same positions along the propagation axis in both runs, and electrons 
are collected from the same plasma volume. Despite differences in minor details, both 
codes consistently reproduce physical details of the self-injection process over the en- 
tire dephasing length. This level of agreement between very different numerical models 
indicates that the results are largely free of numerical artefacts. Importantly, the dis- 
crepancies emerge when the interaction develops noticeable non-cylindrically symmetric 
features, and hence the reduced field description of calder-CIRC loses precision. We 
believe that the agreement between the models may be improved in a straightforward 
fashion (viz. using a larger number of poloidal modes) without significantly reducing 
computational efficiency. 

3.3. Effects of numerical dispersion control 

As described above, simulating upcoming experiments will require economising on com- 
putational cost as much as possible without sacrificing physical accuracy. One means of 
reducing longitudinal resolution requirements, and hence allowing a larger time step, is 
to minimise numerical dispersion through a modified electromagnetic update. Here we 
show how numerical dispersion quantitatively affects the injected electron bunch. 

The immediate effect of numerical dispersion is an unphysically low group velocity of 
the laser pulse. While this effect is difficult to observe directly in the laser pulse because 
of more significant changes in the pulse shape, it can be seen in the electron phase space, 
which is of experimental importance. We examine the initially-injected electrons at the 
point where they have rotated in phase space such that the beam has achieved minimal 
energy spread. The minimal energy spread condition is characterised by the phase space 
of the bunch being roughly longitudinally symmetric and in the shape of a 'IT. We find 
from the perfect dispersion simulation that this occurs after the laser has propagated 
approximately 1.8 mm into the plasma. We show longitudinal momentum spectra and 
phase space at this point for both perfect dispersion and normal dispersion in figure |5j 
We find that with the normal dispersion algorithm, the beam achieves lower energy and 
exhibits higher energy spread than with the perfect dispersion algorithm. We also find 
that phase space rotation has occurred more quickly. 

We also compare the two dispersion algorithms at points of minimal energy spread. 
Without dispersion control, the more rapid phase space rotation causes the minimal 
energy spread to occur after just 1.6 mm of propagation rather than 1.8 mm. We show 
the two phase space plots in figure [6] This comparison is relevant since for applications, as 
one would want to design the system such that the injected beam exits the plasma at this 
point of minimum energy spread ( |Hafz et "a77||201 1| ) . It is clear from these plots that with 
the normal dispersion algorithm, the beam has reached lower mean energy (390 MeV) 
at the minimum energy spread point than with perfect dispersion, where the beam has 
mean energy of 460 MeV. In addition, the normal dispersion case exhibits slightly higher 
energy spread and total charge in the bunch. 

We also compare the longitudinal momentum spectra and phase space at the points 
compared with calder-CIRC simulations in the previous section, namely 960 um, 1.24 mm, 
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Figure 5. Top: Electron longitudinal momentum spectra after 1.8 mm propagation for perfect 
and normal dispersion. Bottom: Longitudinal phase space for perfect dispersion (left) and normal 
dispersion (right). 
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Figure 6. Longitudinal phase space for perfect dispersion (left) and normal dispersion (right) 
at the minimal energy spread point for both algorithms. For perfect dispersion, this is after 
1.8 mm of propagation, and for normal dispersion after 1.6 mm. 

and 2.4 mm. These comparisons are shown in figure [7| We can see, especially in the later 
two plots, that the injected bunch in the normal dispersion simulation shows both lower 
mean energy and greater phase space rotation than in the perfect dispersion run, which 
agrees better with the C alder- CIRC simulations as seen in the previous section. 
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Figure 7. Comparisons of longitudinal momentum and phase space. Momentum spectra for 
perfect and normal dispersion (left), the perfect-dispersion phase space (center), and the nor- 
mal-dispersion phase space (right). The top row shows the electron distribution at 960 um of 
propagation, the middle row at 1.24 mm, and the bottom row at 2.4 mm. 



While these discrepancies are small, they are noticeable and consistent with numerical 
group velocity error. As LPA system designs are refined, and diagnostics and control over 
the laser pulse and plasma improve, it will be important to control numerical effects on 
this level to optimise parameters through simulation. The perfect dispersion algorithm 
allows us to do so while still using low longitudinal resolution for computational efficiency. 



4. Conclusions 

In this paper we have demonstrated the utility of using computationally efficient, fully 
explicit 3D PIC codes to describe and explain the physical phenomena accompanying 
electron acceleration until dephasing in a self-guided LPA in the blowout regime. Electron 
self-injection and its relation to nonlinear dynamical processes involving the laser pulse 
and bubble were explored. Two approaches to reducing the computational cost of the 
simulations were considered. 

First, using the Cartesian code VORPAL with a newly developed perfect dispersion 



algorithm (Cowan et al. 2012), vorpal-pd, made it possible to use large grid spacings 
(~15 grid points per wavelength in the direction of propagation) and proportionally larger 
time steps. This approach reproduces the correct group velocity of a broad-bandwidth 
laser pulse. The red-shift, self-compression, and depletion of the laser pulse were thus 
described correctly, with proper resolution of all important physical scales. 

Second, the well-preserved axial symmetry of the problem allowed us to use a re- 
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duced geometry description, with poloidal-mode decomposition of currents and fields. 
This approach was realised in the code calder-circ ( jLifschitz et a/.||2009| ). By using 
only two modes, we approached the performance of a 2D code, at the same time preserv- 
ing the correct cylindrical geometry of the interaction. The high computational efficiency 
of calder-circ allowed us to use a very high longitudinal resolution (~50 grid points 
per laser wavelength in the direction of propagation) and a large number of macropar- 
ticles (~50 per cylindrical cell), eliminating numerical dispersion and strongly reducing 
the sampling noise. This high resolution simulation did not indicate any new physical 
effects relative to the vorpal-pd runs, and did not exhibit significant differences in the 
quantitative results. Even with strong violation of cylindrical symmetry (such as near 
the dephasing limit, when the pulse was transformed into a two-cycle relativistic piston), 
the CALDER-CIRC results remained qualitatively correct. 

Both codes described precisely the self-focusing of the laser pulse, the oscillations of its 
spot size, and related oscillations of the bubble; electron self-injection into the oscillating 
bubble and formation of a quasi-monoenergetic bunch; laser pulse frequency broadening 
and self-compression into the relativistic piston; constant elongation of the bubble dur- 
ing the piston formation; and uninterrupted electron injection eventually overloading the 
bubble. The codes showed excellent agreement on the locations of initiation and extinc- 
tion of injection, on the collection volume, and on parameters of the quasi-monoenergetic 
component in the electron spectrum, indicating that the results are free of numerical arte- 
facts. It is especially interesting that the calder-circ simulation with just two poloidal 
modes did not lose accuracy and preserved the correct group velocity (agreeing with the 
VORPAL-PD run) even when the laser pulse was compressed down two cycles. 

We thus conclude that (1) using perfect dispersion, taking a coarser grid and larger 
time steps, and using higher-order splines for macroparticle shapes to suppress the sam- 
pling noise, or (2) neglecting high-order non-axisymmetric field and current components, 
thus reducing the dimensionality of problem are both effective and promising means to 
increase the computational efficiency without sacrificing fidelity. Both of these meth- 
ods are applicable to the design of upcoming experiments on GeV-scale acceleration of 
electrons with 100 TW-scale lasers. 
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